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The nature of randomness in disordered packings of frictional and frictionless spheres is inves- 
tigated using theory and simulations of identical spherical grains. The entropy of the packings is 
defined through the force and volume ensemble of jammed matter and shown difficult to calculate 
analytically. A mesoscopic ensemble of isostatic states is then utilized in an effort to predict the 
entropy through the definition of a volume function dependant on the coordination number. Equa- 
tions of state are obtained relating entropy, volume fraction and compactivity characterizing the 
different states of jammed matter, and elucidating the phase diagram for jammed granular matter. 
Analytical calculations are compared to numerical simulations using volume fluctuation analysis and 
graph theoretical methods, with reasonable agreement. The entropy of the jammed system reveals 
that the random loose packings are more disordered than random close packings, allowing for an 
unambiguous interpretation of both limits. Ensemble calculations show that the entropy vanishes at 
random close packing (RCP), while numerical simulations show that a finite entropy remains in the 
microscopic states at RCP. The notion of a negative compactivity, that explores states with volume 
fractions below those achievable by existing simulation protocols, is also explored, expanding the 
equations of state. The mesoscopic theory well reproduces the simulations results in shape, though 
a difference in magnitude implies that the entire entropy of the packing may not be captured by 
the herein presented methods. We discuss possible extensions to the present mesoscopic approach 
describing packings from RLP to RCP to the ordered branch of the equation of state in an effort to 
understand the entropy of jammed matter in the full range of densities from RLP to FCC. 



I. INTRODUCTION 

Granular materials fall under the scope of athermal 
systems, which includes glasses, colloids and gels, among 
others. These athermal systems exhibit non-equilibrium 
behavior, such that equilibrium statistics is insufficient in 
its attempt to describe the system dynamics. These sys- 
tems are thereby considered "complex", and their char- 
acterization finds application in fields from chemistry to 
fluid mechanics and beyond. For granular systems, in 
particular, a phase transition P, [J occurs when gran- 
ular materials are compressed such that they develop 
a nonzero stress in response to a strain deformation 
d, 0, H, H, This transition, referred to as the jam- 
ming transition, occurs at a critical volume fraction, (j)c, 
depending on interparticle friction and preparation pro- 
tocol. Analysis of the jamming transition produces a 
phase diagram of jammed granular matter for identical 
spheres, characterized by 0c and the average mechanical 
coordination number Q. 

The existence of boundaries in the phase diagram of 
Q are related to well-defined upper and lower limits in 
the density of disordered packings; random clo se p acking 
(RCP) and random loose packing (RLP) M- How 
to properly define RCP and RLP remains a longstand- 
ing open question in the field. It has been suggested [ij 
that treating a jammed system via the volume (V) en- 
semble introduces an analogue to temperature in equilib- 
rium systems. This analogue, "compactivity", is a mea- 
sure of how compact a system could be. Within this 
framework P, Q, RCP is achieved in the limit of min- 
imal compactivity and RLP is achieved in the limit of 
maximal compactivity. Therefore, the boundaries of a 



phase diagram for jammed matter could be defined by 
the limits of zero and infinite compactivities. 

In order to approach jammed systems with a statistical 
ensemble approach, a definition of RCP and RLP requires 
proper definitions of jammed states and the concept of 
randomness [ll| . In an attempt to rigorously define 
jammed states, Torquato and coworkers have proposed 
three categories of jamming jT^]: locally, collectively and 
strictly jammed. This definition is based purely on geo- 
metrical considerations and therefore it is only sufficient 
for frictionless grains. Frictional systems incorporate ge- 
ometrical constraints but are dominated by inter-particle 
normal and tangential contact forces [l^. In Fig. [T]we 
see a hard sphere system is not locally jammed if only 
normal forces are considered, since the ball can freely 
move in the vertical direction. The same geometrical con- 
figuration is locally jammed if friction is allowed between 
the particles, revealing the importance of forces in the 
definition of jamming for frictional particles. Therefore, 
a definition of the jammed state for granular materials 
considering only geometrical constraint is insufficient to 
describe frictional grains. 

Frictional systems further exhibit an inherent path de- 
pendency, as granular contacts between grains result in 
the loss of energy conservation. Approaches based on the 
potential energy landscape thereby cannot be used 
for granular materials, as such a potential does not exist 
for the non-conservative frictional contact force. In this 
study, our framework is based on statistical mechanics 
[T^ . defining the jammed state at the V-ensemble supple- 
mented by force and torque balance conditions, wherein 
volume replaces energy as the conservative quantity for 
a statistical ensemble. The free volume associated with 
each particle in the packing is calculated as a function 
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FIG. 1: (a) The inset shows a ball in 2d under mechanical 
equilibrium by two nearest neighbor contacts. The ball is 
not jammed under a normal force interaction. It jams when 
tangential forces are present. 

of the geometrical coordination number using a coarse- 
grained mesoscopic theory of quasi-particles. This allows 
one to define the RLP and RLP states through the en- 
semble of isostatic states. 

Randomness in statistical systems is typically charac- 
terized by the entropy, the equation of state derived from 
the number of microstates available to the system. In 
equilibrium statistical mechanics, entropy provides the 
link between these microstates and the macroscopic ther- 
modynamic properties of the system. We explore that the 
concept of randomness is well-defined for the V-ensemble 
following the Gibbs distribution [14], and is different from 
the measurement of randomness of single packing in term 
of the ensemble of order parameters proposed in 
Therefore, calculating the entropy within the V-ensemble 
can relate the available microscopic volume for each grain 
to the macroscopic system properties, such as volume 
fraction, average coordination number, and compactiv- 
ity, in the case of frictional hard spheres. 

We first investigate frictional packings of equal sized 
spheres at the jamming transition generated via com- 
puter simulations. As the volume fraction approaches 
the jamming transition from above, (f) — > <^^, the sys- 
tem approaches the isostatic condition and observables 
are shown to scale with the distance from the jamming 
transition as a power law oi 4> — 4>c [3, 0, 0|; including 
stress, average coordination number and elastic moduli. 
We therefore consider the jamming transition as the limit 
at which the stress tends to zero and the average coor- 
dination number tends to a finite, non-zero, value. Me- 
chanical equilibrium imposes an average mechanical co- 



ordination number, Z ^ larger or equal than the minimum 
isostatic coordination as co nject ured by Alexander [l^ 
(see also H i, 0, [H, [H, [Hill, [H). The ensemble of 
packings at the jamming transition explores the phase 
diagram of jammed matter by assuming the system is 
exactly isostatic at the transition, and that the isostatic 
condition varies as a function of the inter-particle friction 
coefficient [8|]. 

We compute the equations of state, entropy and com- 
pactivity, as a function of volume fraction, ranging from 
RLP to RCP. The entropy is calculated by two methods. 
First, a direct analysis of volume fluctuations via Ein- 
stein Fluctuation theory explores the clustering of micro- 
scopic volumes. Second, graph theoretical methods using 
the Shannon entropy analyze the network forming prop- 
erties of the granular system. These simulations reveal 
that random loose packings have a higher disorder than 
random close packings. Further, packings approaching 
RLP have a higher compactivity than those approaching 
RCP. 

Then we perform theoretical calculations under the 
quasi-particle mesoscopic approximation of |^], where a 
coarse-graining over a mesoscopic length scale of several 
particle diameters is implied, giving rise to a mesoscopic 
configurational entropy, achieving a minimal value at the 
volume fraction of RCP and maximal value at the RLP 
limit. The results define RCP and RLP at the mesoscopic 
level in general agreement with simulations, suggesting 
that the concept of randomness in [l4| together with the 
notion of the jammed state in Q are useful. The numer- 
ical results further suggest that the mesoscopic entropy 
requires augmentation to include the entropy of the mi- 
croscopic states neglected at the mesoscopic level. 

The maximal volume fraction for jammed spheres in 3d 
created using purely random protocols occurs at (/)rcp — 
0.64. Packings above RCP exist with some degree of or- 
der, up through the perfectly ordered FCC state with 
'/'FCC = 0.74 in 3d. It is of interest to understand how 
one would expand the existing mesoscopic theory to in- 
clude all packings from RLP to FCC, and what effect 
this would have at the transition between disordered and 
ordered packings at RCP. While entropy tends to zero as 
we approach FCC, the existing mesoscopic theory con- 
siders entropy minimal at RCP, accounting only for dis- 
ordered states. It remains an open topic as to whether a 
phase transition occurs at RCP, noted by a discontinu- 
ous change in the equations of state, or whether disorder 
decays smoothly when approaching FCC. We discuss pro- 
pose plausible scenarios to rationalize the transition be- 
tween RCP and FCC as the volume fraction is increased 
by partial crystallization. We speculate that at RCP a 
thermodynamic transition, either continuous or discon- 
tinuous, may occur. Such a transition can be described 
by a full theory that includes both ordered and disordered 
states and is beyond the scope of the present work. We 
stress that this is a hard sphere transition different from 
the jamming transition obtained for deformable particles 
as the external pressure approaches zero. In this work, 
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hard sphere packings are numerically realized by simulat- 
ing soft particles in the limit of zero pressure using the 
"spht" algorithm explained in Section HH A 

Existing packing protocols exploring jammed packings 
may not probe the entire phase diagram for jammed mat- 
ter, as packings with a negative compactivity may exist 
with volume fractions beneath the minimum value [20| . 
The mesoscopic theory is analyzed in an effort to char- 
acterize these packings that are inaccessible via random 
generation protocols. 

II. SIMULATIONS AND RESULTS 
A. Packing Preparation 

First, we investigate the entropy of jammed granu- 
lar matter by analyzing computer generated packings 
of 10,000 spherical equal-size particles (a reader famil- 
iar with Q can refer to section II. C). 

Two spherical grains in contact at positions fi and r2 
and with radius R interact with a Hertz normal repulsive 
force [m and Mindlin tangential contact forces . 

The Hertz force is defined as: 

F„ = ^ KR'^'S'^', (1) 
and an incremental Mindlin tangential force is defined as: 

AFt = kt{R5)^/^As, (2) 

Here the normal overlap is (5 = (l/2)[2i?— i;2|] > 0. 
The normal force acts only in compression, Fn ~ when 
6 < 0. The variable s is defined such that the relative 
shear displacement between the two grain centers is 2s. 
The prefactors fc„ = 4G/(1 - ly) and h = 8G/(2 - ly) 
are defined in terms of the shear modulus G and the 
Poisson's ratio v of the material from which the grains 
are made. We use G = 29 GPa and i' = 0.2 typical 
values for spherical glass beads and we use i? = 5 x 
10^^ m and the density of the particles, p = 2 x 10'^ 
kg/m"^. Viscous dissipative forces are added at the global 
level affecting the total velocity of each particle through 
a term —jx in the equation of motion, where 7 is the 
damping coefficient related to the viscosity of the medium 
T] = 7/(67ri?). These dissipative forces ensure that the 
granular system cannot 'rattle' forever. We measure the 
time in units of to = R^Jp/G, the compression rate in 
units of Fq = b.^t^^ and the viscosity in units of 770 = 
8.2R^p/to. The dynamics follows integration of Newton's 
equations. 

Sliding friction is also considered: 

Ft < fiF„. (3) 

That is, when Ft exceeds the Coulomb threshold, /iF„, 
the grains slide and Ft = fiFn, where fj, is the static 
friction coefhcient between the spheres. 



The critical volume fraction at the jamming transition, 
, is achieved by the "split" algorithm as explained in 
, allowing one to obtain packings at the critical density 
of jamming with arbitrary precision. 

Initially, a dilute particle configuration is generated 
randomly, usually with a volume fraction 0.30 ^ 0.36. 
Then, an extremely slow isotropic compression, without 
friction, is applied to this configuration until the system 
reaches 0^ in an unjammed state. For frictional packings 
the critical volume fraction, (/>c, is such that (pc < 0.64. 
Additionally, the mechanical coordination number, Z, or 
the number of contacts for a given particle which apply a 
force to maintain the jamming condition at (/)c, is Z < 6 
for frictional systems. Therefore, the system now at (pi is 
allowed to relax while maintaining the frictionless condi- 
tion, such that the system is unable to achieve jamming. 
Z and the pressure decay to zero as we relax the system 
below jamming. 

After obtaining the relaxed, unjammed and friction- 
less state with initial volume fraction a particular fi 
is given to the particles and compression is applied with 
a compression rate F until a given volume fraction 
Then the compression is stopped and the system is al- 
lowed to relax to mechanical equilibrium by following 
Newton's equations without further compression. The 
split algorithm searches for cjjc by setting upper and lower 
boundaries for (jjc and dividing the size of those bound- 
aries in half by iteratively compression (or expanding) 
the system, followed by relaxing and then testing for a 
non-zero stress in the system, as outlined in ^SJ. This 
numerical process is repeated for packings with varying 
along the range of /i between and 00. The results 
fill a phase diagram for jammed identical spheres right 
at the jamming transition as obtained in 8] and shown 
in Fig. [3 

The friction coefficient ranges from to cx3 producing 
packings with coordination number varying from Z fa 6 
to Z « 4, respectively. Wc find that there exits a com- 
mon function Z{fj,) over the different F and (pi (see [8]). 
For /i —f 00, (p ranges from the RLP limit 0rlp ~ 0.55 
obtained when F — > and (pi < 0.55 to the RCP limit 
(p^CP ~ 0.64 obtained for larger F and (j)i — > 0.64. 
For fj, ~ 0, the density is approximately (p ~ 0rcp 
with Z ~ 6. All of the packings used herein, along 
with the computer codes necessary to generate the pack- 
ings and calculate their entropy can be downloaded at 
http:/ /j amlab.org. 



B. Phase Diagram 

Simple counting arguments, neglecting correlations be- 
tween nearest neighbors, consider that a necessary con- 
dition for mechanical equilibrium is that the number of 
independent force variables must be larger or equal than 
the number of linear independent force/torque balance 
equations. Alexander [lB| conjectured that at the transi- 
tion point for frictionless spherical packings the 
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system is exactly isostatic with a minimal coordination, 
Z — 2d = 6 in 3d. Such a conjecture can be extended 
to the infinite friction case, where Z = (i+l = 4[ll,0]. 
In the presence of finite intcr-particlc friction coefficient 
/i, there exists a dependency of Z and fi suggested by 
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simulations 

Figure [5] shows the phase diagram used for all equa- 
tion of state calculations presented herein, as obtained 
using the "split" algorithm as described above in the 
plane ((j)^, Z) (for simplicity, in what follows we denote 
4> = 4>c- That is, we understand that all packings consid- 
ered herein are at the jamming transition and are hard 
sphere packings). As discussed in [sj, packings along the 
RCP line for finite /i are most difficult to obtain, most 
notably near C point in Fig[2l resulting in higher values 
of the lowest achievable stress for those packings. The 
G-linc, a.t Z — 4.0, indicates the theoretical Z for infi- 
nite friction packings. The grey line at Z = 4.2 indicates 
the approximate lowest achievable Z possible using the 
present "split algorithm". The solid color lines in Fig. [2] 
are averages used in all following calculations. Near the 
RLP and RCP lines for a some fixed ^, we observe no- 
tably higher values of Z. These values are not included in 
the average. This allows us to use a constant value as an 
approximation for the mechanical coordination number. 

We compare this result to the phase diagram as pre- 
dicted by theoretical model asserted in Q. Note that 
the isostatic condition [l^ predicts Z = 6, while Fig. [2] 
includes packings with 6 < Z < 6.2, and (f) > 0.634 as 
predicted by the theory for RCP. We suggest that these 
packings are new microstates of jammed matter (indi- 
cated by the shaded portion of the phase diagram) which 
are not accounted for in the mean- field version of [1, [2^ . 
While they remain a component of the ensemble gener- 
ated via the above described simulation protocol, their 
existence remains a topic of ongoing study. 
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FIG. 2: The phase diagram of jamming from simulation re- 
sults. We use the same algorithm to generate the packings as 
done in [|. Here, we use 10,000 particles, while in Q only 
1000 were used. A larger number of particles is necessary 
for accurate entropy calculations. The volume fraction is de- 
noted by as opposed to for simplicity. Horizontal lines 
show the average coordination number used for packings of 
constant /x. The dashed line represents the theoretical RLP 
line. The solid vertical line ai tf) — 0.634 is the theoretical 
RCP line obtained in Notice that some packings exist to 
the right of the RCP line. Such packings are not captured by 
the theory, indicating that microscopic fluctuations beyond 
the mesoscopic theory of @j are important close to the RCP 
state. The solid grey line at Z = 4.2 indicates the lower limit 
for Z available using the present split algorithm. The J-point, 
located at {(f>,Z) — (0.634,6), is the theoretical frictionless 
jamming point. The L-point, located at {4>,Z) = (0.536,4), 
is the theoretical jamming point for /j. ^ 00 with A — > 00. 
The C-point, located at {(p, Z) = (0.634,4.0), is the theoretical 
jamming point for /x — > 00 with A — > 0. The G-line, Z = 4.0 
is the theoretical average Z achieved for all infinite friction 
packings of identical spherical grains in 3d. 



C. Entropy from Voronoi Volume Fluctuations 

In the absence of energy conservation, a different sta- 
tistical approach is necessary to describe the ensemble 
properties of jammed granular matter. Along this line of 
research, Edwards [14] proposes replacing the system en- 
ergy by the volume as the conservative quantity such that 
a microcanonical partition function of jammed states can 
be defined and a statistical mechanical analysis is plausi- 
ble. Therefore, a microscopic volume must be associated 
with each grain. 

As detailed in Jamming I [2^, the definition of a 
Voronoi cell is a convex polygon whose interior consists 
of all points closer to a given particle than to any other. 
Further, it is additive and tiles the system volume com- 
pletely. The formula for the Voronoi volume of a particle, 
i, in terms of particle positions for monodisperse spheri- 



cal packings in 3d is f23| 

wvor ^ i f J_ ds, (4) 

3 J \2R 3 cos9.,J ' ^ ' 

where rij is the vector from the position of particle i to 
that of particle j, the integrand is over all the directions 
s forming an angle 9ij with r^, and R is the radius of 
the grain. The Voronoi volume is used to tile the to- 
tal system volume, and replaces energy as the conserved 
quantity in a new micro-canonical ensemble for jammed 
granular matter. Therefore, fluctuations in Voronoi cell 
volumes are related to the compactivity of the jammed 
system, much like energy fluctuations are directly related 
to the system temperature in equilibrium thermodynam- 
ics. We notice that the Voronoi-Delaunay decomposition 
is the basis for Hales proof of the Kepler conjecture [2^ . 
Below, we treat the monodisperse case. Other cases will 
be treated in subsequent papers. 

Next, we calculate the entropy of the numerical pack- 



ings in Fig. [2] from Voronoi volume fluctuation analo- 
gous to Einstein Fluctuation theory. We first define the 
Voronoi cell associated with each particle i and calcu- 
late its Voronoi volume Wj. Calculation of a Voronoi 
cell volume begins by defining the polygon between two 
Delaunay contacts, having finite number, m, vertices. 
Two grains are considered Delaunay contacts if their cor- 
responding Voronoi cells share a face. Delaunay con- 
tacts are determined by the network of grain positions 
and radii calculated using QHuU software, available at 
://www. qhull.org. The contribution of this polygon 
to calculating the Voronoi volume comes from the abil- 
ity to associate a pyramid, comprised of the center of 
each particle as the apex, and the m-sided polygon as its 
base, as shown in Fig. [3] (schematically in 2d for simplic- 
ity), to each particle. The two pyramids are symmetric. 
The volume of this pyramid is the contribution to the 
Voronoi volume of the cell surrounding a particle, ex- 
clusive to the particular Delaunay contact which shares 
the polygon base. Repeating this process for each De- 
launay contact results in the complete Voronoi volume 
surrounding a particle. The Voronoi volume is thereby 
the microscopic volume associated with each grain. 

We perform statistical analysis of the volume fluctua- 
tions by considering a cluster of n particles. The Einstein 
fluctuation relation is defined as follows [1^ [2^ : 

al = {(Wn-{Wn))^)^XX^d{Wn)/dX (5) 

Equation ^ is analogous to equilibrium thermody- 
namics, replacing energy and temperature by volume and 
compactivity, X, in the Edwards picture. Note that A is 
the analogue of the Boltzmann constant ks that defines 
the units of compactivity. 

We calculate the average volume, (Wn) and fluctua- 
tions CT„ = {(Wn — (Wn))^), where (•) is an average over 
many n-clusters. From the large n behavior we extract 
the fluctuations versus volume fraction, for every pack- 
ing depicted in FigO Figure S] shows the fluctuations as 
a function of n for packings with inflnite friction, display- 
ing the largest range of volume fractions in the data used 
herein. We find that for sufficiently large n ^ ric, the 
fiuctuations scale with n and therefore are extensive and 
well-defined. 

Figure [5] shows the approximate value of Uc at which 
the extensive nature of the fluctuations reaches its max- 
imal value, as a function of (j). Lower volume fractions, 
approaching RLP, require higher values of n w 1000 to 
achieve this condition. This result contrasts with the 
results of [27], although there the system was smaller, 
N w 100, and two-dimensional. Reference [l^ acknowl- 
edges that if grain volumes can be treated as indepen- 
dent random variables, then the fluctuation in clusters of 
n Voronoi cell volumes should scale with n. For clusters 
of jammed grains, this was not observed in [27| . indicat- 
ing the existence of correlations between the Voronoi cell 
volumes within a cluster. However, in the present study, 
the value of n at which the fluctuations are extensive is 




FIG. 3: (a) Example of 2d Voronoi cell - All of the calcula- 
tions are done in 3d but are shown here in 2d for simplicity. 
The line between the centers of particles i and j is defined by 
ij, equivalent to r^j. The line perpendicular to the bisection 
of ij is defined by ab, intersecting ij at point c. 4 additional 
particles are Delaunay contacts of particle i, such that a 5 
sided polygon (pentagon) surrounds particle i, defining the 
Voronoi cell of particle i by virtue of intersecting bisecting 
lines between each pair of Delaunay contacts. Points a and 6 
are defined as the boundary of the Voronoi cell line between 
particles i and j. A triangle is thereby formed by points i and 
ab, the area of which is the contribution of the Voronoi cell 
of particle i exclusive with its Delaunay contact to particle 
j. This process is repeated for all Delaunay contacts of i to 
give the entire Voronoi cell area. Note that a symmetric area 
to Aiab exists as Ajab, and can be appfied to the Voronoi 
cell area of particle j. (b) Example of 3d Voronoi cell - The 
line between the centers of particles i and j is defined by ij, 
equivalent to rij. The plane perpendicular to the bisection 
of ij is defined by p, intersecting ij at point c. Note that 
particles i and j are identical spheres, with particle j appear- 
ing smaller only to illustrate the 3d properties of the system. 
Plane p is intersected by m other planes, creating an m-sided 
polygon between particles i and j. Each plane intersecting 
plane p (not shown) is a plane bisecting ik, the line between 
the centers of particle i (or j) and another particle k in the 
system, where k is one of m particular particles. A pyramid 
is thereby formed using the m-sided polygon as the base, and 
i (or j) as the apex. This pyramid is symmetric over plane p, 
and its volume is the contribution to the Voronoi volume of 
particle i from particle j, or vice versa, exclusively. The vol- 
ume of the pyramid is calculated by separating the pyramid 
into 8 smaller pyramids, using the triangle composed of one 
of the m available sides, and c as its base and i as its apex. 
This is illustrated by using ab and c as the base of a pyramid 
with apex i. The volume of this pyramid is calculated and 
the process is repeated for each of the m sides, adding each 
obtained volume to the Voronoi volume of both i and j. The 
entire process is then repeated for all Delaunay contacts for a 
given particle, resulting in the total Voronoi volume for that 
particle. 



a function of (j), shown in Fig. [5l In [23] , this phenomena 
is observed, but the density is apparently independent of 
the volume fraction, and occurs at the same value of n for 
each packing. For packings approaching RCP, the exten- 
sive nature of the fluctuations occurs at ric ~ 100, lower 
than Tic ~ 1000 for RLP. Figure [5] shows the fluctuation 
density, ((W„ - (W„))2)„ = ((w„-(w„))^) ^ fluctuation 
per grain, as a function of (f> for all packings used herein 
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FIG. 4: {AWl) versus Uc for packings with ^ — > cx3, with FIG. 5: Wc vs (j) showing maximal value of extensive nature. 
Z « 4.3, for different (j). Values of 4> are taken along fj, — > oo, with Z ~ 4.3, to display 

largest range of (j). 



obtained for n > ric- This data represents the main equa- 
tion of state for all of the numerical packings considered 
herein. Each color curve represents packings with a fixed 
Z (or /i) as indicated in the figure. We note that while 
the fluctuations for all Z{ii) in this study collapse onto 
a single curve, as illustrated in Fig [6l the limit of inte- 
gration, </)RLp(-Z') in Eq. Q, changes as discussed in the 
phase diagram of [8] , increasing as fj, decreases and effec- 
tively depending on Z as 0rlp(^)- Indeed, (j)Rhp{Z) is 
the theoretical RLP line depicted as a dashed line in Fig. 
m These equations of state should be compared with an 
analogous equation of state obtained in j26| for jammed 
packings equilibrated using a fluidized bed. The fluc- 
tuations presented in Fig. [6] monotonically decrease as 
(j) increases, and do not display a parabolic shape as de- 
picted in [2^]. While the protocol used herein, the "split" 
algorithm presented in Section |TT1 A and the protocol 
for the experiments of [1^ differ, we would expect that 
the equation of state should be the same. Elucidation 
regarding this difference requires further investigation. 

An important note is that this extensive relationship 
occurs well before n is large enough such that finite size 
effects of the system force the fiuctuations to tend to zero. 
Further, the linear relationship extracted from n-clusters 
is different from that extracted from n randomly chosen 
Voronoi cells, implying a correlation between Voronoi cell 
volumes, revealed using the clustering technique. Analy- 
sis of the fluctuation density reveals the following formula 



al + (W„AW„+i) + (W„_iAW„) 
(W„)(AW„+i) - (W„_i)(AW„), 



(6) 



where cr^ is the single particle, or microscopic, fluctua- 
tion in Voronoi volume, and AWn+i = W^l\ - {W""'') 
is the n -t- 1 free Voronoi volume to be added to the clus- 
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FIG. 6: {AW^)„ versus (j) for packings with different fric- 
tion coefficients producing different mechanical coordination 
number versus the volume fraction. Each dataset with fixed 
coordination number corresponds to the packings along each 
horizontal line in Fig. [2] 



ter of n free Voronoi volumes being averaged. If volumes 
are chosen randomly, as opposed to the clustering condi- 
tion, the fluctuation density collapses to the microscopic 
fluctuation density, af, as the correlation tends to zero. 
Therefore, the fiuctuations will scale exactly with n, as 
indicated in [S^. Further, if the averaging process is 
taken to be equal to or larger than the system size, then 
(AW„+i) = {AWn) = 0, such that Eq. ® is rewritten 
as 



Aal ='^1 + {WnAWn+l) + {Wn^lAWn) 



(7) 



Equation ([7]) thereby provides an analytical form for 
the curves presented in Fig. 31 as a function of n. Future 
work may take into consideration nth neighbor coordina- 
tion and distance distribution, as in 1281. 
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Compactivity. — The compactivity is then obtained via 
the integration of Eq. ([5]) : 



0RLP (Z) 



d{Wn 



(8) 



where we use that (j){X oo) — > 0rlp [Sl- Since Voronoi 
volumes are additive, (yV„) = (W) = NVg/(j). The fluc- 
tuations in Voronoi volume are divided by the number 
of grains, N , thereby introducing the fluctuation density, 
shown in Fig. [S]into Eq. ([5]). Therefore, the above inte- 
gration is rewritten as: 



4>(X) 



d(j) 



(9) 



and we may then utilize the fluctuations as a function of 
(/>, and integrate along a line of constant Z{fi). The as- 
sumption that a'i(V)/(V)'^ = (t^(0)/((/))^ is not utilized 
here, as done in [20|, explaining the different functional 
form for the compactivity of Eq. ^ from that of [l^ . 

Following the above presented method, (/)rlp (Z) is ex- 
tracted from the phase diagram, and used as a limit of 
integration in order to calculate X{(j)) from fluctuations 
in Voronoi volume, in Eq. ([5]). This introduces the de- 
pendence on Z as (t>{X) is obtained by integrating Eq. 
Q numerically by applying a fitting function to the nu- 
merical data of Fig. [SI The equation of state, (l){X) for a 
given Z, is plotted in Fig. [7]for different values of the av- 
erage coordination number of the packings, Z(fi), reveal- 
ing that as we approach 0rcp ~ 0.645, X ^ 0, regardless 
of the value of /i. Further, X — s- oo as we approach ((irlp, 
with the smallest volume fraction of the RLP appearing 
for /i — > cx) and Z w 4 in the high-compactivity limit, 
0RLP ~ 0.55. The compactivity curve plotted in Fig. [7| 
is continuous, even though the volume fluctuation data 
is the result of a discrete set of simulations as seen in 
Fig. [6l This is due to the fact that wc smoothly inter- 
polate a continuous curve for the volume fluctuations as 
a function of (f), resulting in a smooth integration for the 
compactivity, and subsequently the entropy through Eq. 
®. 

Entropy. — The entropy, S, and its density, s — S/N, 
are obtained by integrating 



X- 



dS_ 
dV' 



(10) 



By virtue of having a fixed total volume, V, for any 
particular system defined by 0, we can substitute V = 
NVg/(t) such that: 



{XlVg)-' = 



ds 



(11) 



Using the concept that 0rcp is a fixed value in the 
phase diagram for all values of Z(fi), we integrate be- 
tween the limits of (/)rcp and the desired 0. 
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FIG. 7: <j) versus X from the integration of Voronoi volume 
fluctuations. The smoothness of the curves is due to the fact 
that we use fitting functions for the data in Fig. [6]to perform 
the integration of Eq. © 



s{(f>) - s((/'RCp) = A 



(12) 



Equations ([9]) and (fT2|) can be combined to provide 
an equation for s as a function of the Voronoi volume 
fiuct nations. 



s(0) = A 



0RCP 



RLP ( 



5RCPJ 



(13) 

Therefore, we can calculate s(0)/A, based on the fiuc- 
tuations of a packing of jammed grains for a fixed Z, 
following the horizontal lines in Fig. [51 Integration of 
the compactivity curve achieved via simulation provides 
the entropy, up to a constant value at 0rcp, as defined 
by Eq. The entropy of the packings from simula- 

tions in Fig. [21 is plotted in Fig. [5| as a function of (j) 
for different values of Z. The shape of the entropy curve 
is in qualitative agreement with that calculated by Aste 
and collaborators using X-ray tomography techniques to 
determine the position of particles inside large scale pack- 
ings, as shown in (29| . 



D. Entropy from Information Theory 

Analysis of the entropy from fluctuations in Voronoi 
volume clusters provides a value for entropy density up to 
a constant of integration s(0rcp), as shown in Eq. (|13p . 
To obtain the entropy of RCP we use an independent 
method based on information theory as developed in [s^, 
HH, related to the thermodynamic entropy in [s^l, and 
applied to emulsion systems in (ss'l which does not require 
a constant of integration. This method provides a second 
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FIG. 8: Entropy density versus (p from the integration of X. 

estimate of the entropy for all volume fractions to be 
compared with Fluctuation theory. 

We use the Voronoi cell and Delaunay triangulation for 
each particle to define a Voronoi network by considering 
contacts when a Voronoi side is shared between two par- 
ticles, and hence are Delaunay contacts, as shown in Fig 
[3l In order to facilitate periodic boundary conditions, we 
surround the finite box enclosing all Voronoi cells with 
26 virtual boxes. These boxes enclose virtual particles, 
translated in all possible combinations from the real box. 
QHuU calculates all Delaunay contacts, and only those 
pairs of contacts which contain at least one real particle 
are considered, while pairs of virtual contacts are dis- 
carded, thus ensuring complete periodic boundary con- 
ditions. 

A graph is constructed as a cluster of n particles that 
are Delaunay contacts, and by means of graph automor- 
phism [s^l can be transformed into a standard form or 
" class" i of topologically equivalent graphs that are con- 
sidered a state with a probability of occurrence pi . Ref- 
erence |34] provides the code to deal with isomorphic 
graphs, which is essential for counting different graph 
classes. The topological equivalence accounts for graphs 
with varied Voronoi cell sizes and locations while retain- 
ing identical lists of consecutive subsets of Delaunay con- 
tacts. In practice, we determine pi by extracting a large 
number m of clusters of size n from the system and count 
the number of times, fi, a cluster i is observed, such that: 

p^ = f^/m. (14) 

A simple example of graph classes can be seen when 
n = 3. Figure [9] shows exactly 4 possible graphs that can 
be achieved, consisting of 0,1,2 and 3 connections be- 
tween 3 Delaunay contacts. Each graph falls into a class 
of i = ^ 3. The random process of selecting a graph 
with n = 3 is done m times, continuously calculating the 
probability of fi, as defined in Eq. HH). 

The Shannon entropy of a clusters of size n is thereby 
defined as: 



n =3 

A A /. 

=3 i=2 i.1 i.O 



FIG. 9: Classes of graphs for n=3. There are only 4 possible 
graphs, contacts, 1 contact, 2 contacts, or 3 contacts. 



= -A^pjup,. (15) 

Equation psp reduces to the thermodynamic entropy 
if one replaces the probability, pi, with the Boltzmann 
factor. The name Shannon entropy does not imply that 
Eq. (|15[) is a new kind of entropy. It merely implies 
that we will calculate the entropy of the packing using 
information theory methods, which is routinely applied 
to sequences. 

Ideally, m is large enough such that the values of pi 
converge to a fixed value for each value of n, such that 
the Shannon entropy converges as well. Computationally, 
for large n, it is not possible to reach this convergence 
within a reasonable time, and certain approximations can 
be applied as outlined in [3^ to decrease the number of 
iterations necessary, summarized briefly here. 

As one iteratively selects clusters, increasing m by one 
each time, new classes of graph are obtained, such that we 
have a total number of different graphs k, where k < m. 
While Fig. [5] depicts the relatively simple case of n = 3, 
higher values of n have a significantly larger number of 
graph classes. Therefore, Eq. (|14p is only an approxima- 
tion of the true probability of observing a cluster i and 
can be rewritten as 

Pi ~ fi/m, (16) 

where the equation becomes an equality as m ^ oo. We 
must therefore replace Eq. by 

i/*(n) = -AV ^In^. (17) 

Graphs with fi jm smaller than 1 jm will likely be ob- 
served only once, if at all. This finite-size effect grows 
with increasing ri, as graphs become very complex. The 
quantity iJi (n) is defined as the contribution to H* (n) 
of the topologies observed once. Measurements where 
Hi{n) exceeds a threshold percentage of H*{n) are not 
considered to be valid measurements. 

In an effort to improve convergence, and thereby de- 
crease simulation time for the calculation of the config- 
urational entropy, a finite sample correction is applied. 
The details of this correction are presented herein, as 
well as in (ssj . 
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Referring now to H* from above, the probability Pi{fi) 
that a certain state i wih be observed exactly fi times is 
given by the binomial distribution 



P^if^)-(Jy.'{^-Pir'^'■ (18) 

Define the probability of a certain event to be pi, the 
number of observed events to be m and the number of 
uniquely observed events to be k. Furthermore let F(x) 
be a function that can be Taylor expanded around pf. 



F{x) = F{p,)+F'{p,){x~p,) + -F"ip,)ix-p,f + ... (19) 

The binomial distribution is concentrated around the 
average (fi/m) = pi, and we obtain the following useful 
approximation using the definition of the variance in a 
binomial distribution: 



(F(^)) = F{p.) + F'ip.){^ -p.) + \F"{p.){{l^-p,f) 
m m 2 m 



FiP.) + lF"(p.f-l^^ 
2 m 



(20) 



Let F{x) = — a; log a;, an obvious choice considering the 
form of the Shannon entropy. Then, F"(x) = —1/x and 



TO 



-pi log Pi 



1 -P^ 

2to 



(21) 



Therefore, 



{H* 



= Y.{fC-)) = -Y^pdogi 



m 

H k-l 



E..(i-k) 

2to 



+ 



A 2to 



(22) 



The Shannon entropy is then the calculated entropy 
from the average of H* obtained from multiple simula- 
tions, plus the binomial correction term (fc — 1)/2to, since 
the sum over all pi is unity, a term clearly tending to zero 
as TO tends to infinity. This approximation works well un- 
der two conditions. First, mp > 1, in order to allow the 
Taylor expansion to converge. Second the contribution 
of the probabilities below l/m to the Shannon entropy 
must be a practically insignificant term. If this is not 
the case, additional binomial correction terms must be 
used. However, these terms will not be independent of 
Pi, therefore making the calculation significantly more 
complicated. It is of interest in the present work to en- 
sure that both conditions necessary for the use of only 
the first term of the binomial correction are applicable. 



Assuming the above described binomial distribution 
of the probability of having a cluster of class i observed 
exactly fi times, the first order finite sampling correction 
to H* therefore results in 



H{n) H*{n) k~l 



X 



A 



2 TO 



(23) 



The entropy must be further corrected to account for 
the Shannon entropy as measured for a crystal struc- 
ture by using the methods presented herein. The en- 
tropy of a crystal structure where N ^ oo should be 
zero, and as such an FCC packing should have zero 
entropy. However, when applying the graph counting 
method explained above to a finite system, a crystal 
structure will present a non-zero entropy. It is impor- 
tant to subtract the entropy of a finite crystal structure 
from the Shannon calculations. In studies of network 
forming materials outlined by [Slj, an empirical correc- 
tion term of g{n) = {d — 1) log(n) is subtracted from 
each value of H(n), where d is the dimensionality of the 
network. The functional form of g{n) is obtained by ap- 
plying the above described Shannon entropy calculation 
directly to an FCC packing, and its results are shown 
in Fig. [TOl During the process of randomly selecting to 
points in the network, some points will inevitably fall in- 
side of the grain boundaries. Points approaching the cen- 
ter of a grain will have a different set of n nearest grain 
centers then points chosen outside of the grain boundary. 
This discrepancy is accounted for via g{n). In the present 
study, this term is augmented by using the exact values 
obtained for g{n), and not the empirical form, as shown 
in Fig [TOl It should be noted that the empirical form 
for g{n) comes very close to the exact values, differing 
only by an approximately exact constant. For calcula- 
tions of entropy density, differences between successive 
values of H{n) result in the cancellation of this constant, 
exemplified below. 

We redefine the Shannon entropy as follows: 



H'{n) H{n) 



9{n), 



A A 

and the entropy density is 

s= lim [H' {n + 1) - H' {n)] 



(24) 



(25) 



Figure [TT] shows the calculation of H'in) versus n for a 
typical packing with ^ = 10000 and 4) = 0.64. We show 
that s{n) converges so rapidly, as shown in Fig. 1111 such 
that even moderate values of n (n > 8) are enough to ob- 



tain a sufficient approximation of s 31| . The calculation 
is repeated for all the packings and the Shannon entropy 
density is plotted in Fig. [T^] versus cj) for different Z{ij). 

When examining values of the Shannon entropy for 
values of (/) < 4'b.cp, Fig. [12] displays an increase in 
the entropy as volume fraction decreases, similar to that 
of the entropy as calculated via Voronoi volume fluc- 
tuation density in Fig [S] However, this increase does 
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FIG. 10: Entropy of crystal (FCC) packing as calculated us- 
ing graph theoretical methods. While the entropy is approx- 
imately equal to (d — 1) ln(n), minus a meaningless constant, 
exact values from simulations are used herein. 
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H'(n)A for n= 10000, <t) = 0.64 
Linear Fit to obtain s{n) 
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FIG. 11: Shannon entropy function of n, with n = 10000 and 
(j) = 0.64. The red line displays a linear fit between n = 8 and 
n = 14, from which the entropy density is extracted. This 
process is repeated for all packings used herein. 



not appear dependent on the mechanical coordination 
number of the packings. Further, the Shannon entropy 
does not increase by the same magnitude as observed 
in the entropy from fluctuations. For instance, from 
the fluctuation theorem calculations of Fig. [8] we find 
srlp — srcp = 1-7A while from the Shannon entropy 
of Fig [T2l we find srlp — srcp = 0.35 A. As suggested 
in [35|, this discrepancy may be due to the additional 
entropy arising from freedom to move grains within the 
packing without disrupting the Delaunay network, and 
hence not affecting the probabilities in Eq. p^ . Analy- 
sis of such volume contributions to the Shannon entropy 
requires a Monte Carlo simulation, in which the available 
phase space volume that the packings can explore for a 
fixed graph is probed [ssj , with the additional constraint 
that the packing must maintain the Delaunay contact 
network under all possible rearrangements. This calcula- 
tion is considered in ^] and will be the topic of future 
study. 

Figure [13 shows that as we approach (pucp for all val- 
ues of Z using information theory, s/X ~ 1.1. We there- 
fore define s{(t>RCp)/X = 1-1, the value of the entropy as 
calculated via graph theoretical methods. 

Thus, the Shannon entropy density provides an estima- 
tion of the entropy for the RCP state, s((^rcp) ~ I.IA, 
serving as the constant of integration for the entropy den- 
sity as realized by Fluctuations Theorem. Under this ap- 
proximation, we can shift the Fluctuation Theorem en- 
tropy of Fig. [H] vertically by srcp = 1-lA as calculated 
via Shannon entropy methods. Figure [01 shows the en- 
tropy shifted by this constant value, as calculated via 
simulation. It is important to emphasize that s(</>rcp) is 
approximate, due to neglecting fiuctuations of the fixed 
Delaunay network as explained above. Nevertheless, we 
remark that the obtained value of srcp is compatible 
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FIG. 12: Shannon entropy density, s{<j)) for all packings used 
herein. The minimum value of the entropy density is achieved 
at RCP, and is used as a constant of integration for the en- 
tropy obtained from Fluctuation theorem. 



with other estimates which found the entropy of the order 
A (see for instance the calculation of the analogous com- 
plexity, S, by Zamponi and Parisi who found S/A ~ 1 



III. THEORETICAL MODEL 



In this section we develop a theoretical framework to 
understand the numerical results in light of [8] . 
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frictional hard spheres is: 
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FIG. 13: Entropy at RCP achieves a value of 1.1 as calculated 
by the Shannon Entropy at RCP and serves as a constant of 
integration for the entropy from Fluctuation theorem in Fig. 
(d} producing the entropy shown in this figure. 



A. Statistical mechanics of frictional hard spheres 



Experiments of shaken grains, fluidized beds and oscil- 
latory compression of grains [23, '2^, [13, [13, [38^ indicate 
that granular materials show reversible behavior, and the 
analogue of the conserved energy, E, in thermal systems 
is the volume V — NVg/cj), for a system with N grains 
of volume Vg at position f^. Thus, the number of con- 
figurations, ri, and the entropy in the micro-canonical 
ensemble of jammed hard spheres is defined as [l3 |: 



n{v) 



S(v-Win)^ ejam(rl) Vr,. (26) 



Analogous to the temperature in equilibrium system 
dE/dS — T, the "temperature" in granular matter is 
the compactivity, X = dV/dS. Here 6jam(^i) is a con- 
straint function restricting the integral to the ensemble of 
jammed states, yV{fi) is the volume function associated 
with each particle taking the role of the Hamiltonian in 
thermal systems which is defined in terms of the Voronoi 
volume in Eq. ([4]). The crux of the matter is then to 
properly define 0jam and W to calculate the entropy and 
volume in the ensemble of jammed matter. 

Volume ensemble. — A minimum requirement of 
0jam('^i) is to cnsure touching grains, and obedience to 
Newton's force and torque laws. As in the numerical 
simulations, the volume function, W(fi), is taken as the 
volume of the Voronoi cell associated with each particle 
at position , for which the analytical form has been ob- 
tained in Eq. Q 8]. The entropy in the V-ensemble of 



'3' 



(27) 

the normal inter-particle force is 
Fiji, the tangential force: f^J' = \fij — 

All quantities are assumed properly a- 
dimensional for simplicity of notation. The terms inside 
the brackets {•} correspond to the jamming constraint 



where 

f N 
J ij 



U — 



function Gjam in Eq. 



and therefore define the en- 



semble of jammed states. The first three (5— functions in- 
side the big brackets impose Newton's second and third 
law. The Heaviside 0— function imposes the Coulomb 
condition and the last (5— function the touching grain con- 
dition for hard spheres, assuming identical grains of unit 
radius. Integration is over all forces and positions which 
are assumed to be equally probable as in the fiat average 
assumption in the micro-canonical ensemble. 

An extra term should be added as 

(2V ^y ) ~ ^) ■^here, for the isotropic 

case, the stress tensor is cTq^ = pSai3j with p the pres- 
sure. Since we are treating hard spheres, the pressure 
p just controls the mean value of forces and does not 
contribute to the statistics. Thus, this term is not 
needed, which means that the angoricity, A — dp/dS, 
is irrelevant for hard spheres. At the isostatic limit 
p — > 0, and further assuming a Boltzmann distribution of 
pressure [10, [U, [H, \^ , similar to that of mesoscopic 
volumes Q, such a term would tend towards unity in 
a partition function for jammed matter. However, the 
angoricity should be considered in the case of deformable 
grains, a system of future studies. 

Clearly, Eq. ([27]) is almost intractable from an ana- 
lytical point of view. However, under the quasi-particle 
approximation of we can define the configurational 
entropy at the mesoscopic level using a corollary of the 
force-balance ensemble: the isostatic conjecture and a 
coarse-grain volume function in terms of the coordina- 
tion number as we show in Section [HU B. 

It is of interest to determine if the ensemble defined 
by Eq. (1271) satisfy the definition of jamming given by 
Torquato [l2| . A definition of jammed configurations 
based on force/torque balance is a necessary condition 
for mechanical equilibrium of a packing with friction and 
frictionless grains. A force-balanced packing is defined as 
the existence of a set of forces {/ij,V balls i in contact 
with ball a} such that the sum of the forces/torques for 
each particle is zero. 



E4 = 0' E/' 



0, 



(28) 
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with the non trivial contact forces '^ij \fij\ = 1 , fij • nij > 
where riy is the unit vector of the contact j. Other 
necessary mechanical constraints should be added too, 
e.g. fij A Uij = for frictionless packings. 

While the above conditions are necessary for jamming 
they are not sufficient. To define a more restricted force- 
balance condition we consider that the contacts around 
a ball i are not degenerated if 3 a neighbor a such that 

/I • fc 7^ 0, Vfc ^ 0. (29) 

This condition assures that at least one contact force is 
off-plane. A packing is restrict force-balanced if it is force- 
balanced and the contacts around all the balls are not 
degenerated. 

A question is raised whether the restrict force-balance 
condition for the frictionless case fits the geometrical def- 
inition of Torquato . It can be proved that packings 
with ((28|) and ((29|) are at least locally jammed. To see 
this, let us formalized the locally jammed condition as 
follows: for any ball i, there is at least one contact fiia 
such that fiia-k > 0, Vfc 7^ 0. Since the contact force fia is 
parallel with riicc for frictionless packings, the condition is 
equal to fia - k < 0. From the non-degenerated condition, 
there must be a contact satisfying fi^- k^ 0. Assuming 
fia ■ k > (otherwise we prove the result directly) , there 
must exist another contact a' that fia' • fc < because of 
the force-balance condition, then fia ■ k = 0. 

Thus we prove a sufficient condition for local jammed 
configuration from the force balance point of view. Ac- 
cordingly, the non-degenerated condition is necessary for 
the proof above. Since it not necessary it could also sat- 
isfy more stringent geometrical definitions, e.g. collec- 
tive jammed, or even strict jammed, but this has to be 
figured out in further investigation. Although the re- 
strict force-balance condition is not always satisfied in 
our simulations, the bias is very small. In frictionless 
packings it implies that the coordination number must 
be larger than d; a small fraction of particles with Z < d 
are found in simulated packings as well as experimental 
ones. By removing the degenerated balls recurrently (for 
instance it is a common practice to remove fioaters), we 
end up with a packing satisfying the constraint of non- 
degeneracy and therefore being locally jammed. Thus 
we expect that many experimental packings satisfy the 
restrict force-balance condition. The force-balance con- 
dition can be extended to the frictional cases, by sim- 
ply adding the constraint between normal and tangential 
forces. Finally, the ensemble in Eq. (I27p is assumed to 
satisfy the restrict force-balance condition. Thus, it is 
understood that the condition of Eq. ((29|) is implicit in 
the ensemble average of Eq. ([?f|) . 

Assuming that the conditions specified in Eq. ([77]) are 
met in the numerical packings, the simulation results can 
be interpreted as the ensemble average Eq. ([27]). How- 
ever, there is a further important distinction between Eq. 
(P7)l and the numerical calculations. Equation ^7} as- 
sures that all configurations at a given volume have the 



same probability. This is the flat average assumption in 
the micro-canonical ensemble that allows for the devel- 
opment of statistical mechanics. Without this assump- 
tion statistical calculations are impossible to perform (see 
however (4^ for a thorough discussion) . There is no rig- 
orous proof that the flat average assumption is correct 
in equilibrium statistical mechanics. Still, its validity is 
widely accepted. For granular matter, use of the flat aver- 
age is much more controversial. Earlier simulations [38| 
indicate some evidence for ergodicity. The ergodic hy- 
pothesis implies not only the equal probability of states 
but also that for sufficiently long times the phase trajec- 
tory of a closed system passes arbitrarily close to a man- 
ifold defined by a constant volume (or energy) [i^. Ex- 
periments indicate reversible behavior, supporting that 
the fiat average can be applied to granular matter, al- 
though this assertion is certainly not true in general. 

We notice that Eq. ([?f|) is difficult to solve. Analyti- 
cal progress can be done by considering a coarse-graining 
of the Voronoi volume function and working with quasi- 
particle theory developed in Q to obtain a configura- 
tional entropy at the mesoscopic level. 

B. Volume Function 

The mesoscopic theory presented in [1] and [11] coarse 
grains the Voronoi volumes of a jammed granular pack- 
ing over a mesoscopic length scale and calculate an aver- 
age volume function. The coarsening reduces the degrees 
of freedom to one variable, the geometrical coordination 
number of each grain, z, that allows for an analytical so- 
lution of the partition function. We find a mesoscopic 
free volume function Q: 

_ jWr) - Vg 2^/3 
w{z) = = , (30) 

Vg Z 

Here, we note that W™'' has been rigorously defined in 
Eq. ([4]), validated its equivalence to the Voronoi volume, 
and its use in comparison between simulation results and 
a statistical mechanics formulation utilizing the Voronoi 
cell as the microscopic volume associated with each grain. 

If the system is fully random we can extend the as- 
sumption of uniformity from the mesoscopic scale to the 
macroscopic scale, such that we arrive to an equation of 
state relating <^~^ = w + 1 with z as: 

(31) 

z + 2x/3 

C. Partition Function 

Below, we briefly discuss the results of regarding 
the phase diagram. 

We assume that the sum over each quasiparticle with 
volume w(z) is the total volume [2^, then the canonical 
partition function for a single particle can be written as: 
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-2iso = giz)e dz 



(32) 



Equation ([32)1 is the single particle partition function, 
such that the full partition function for N particles is 

The density of states, g{z) is assumed to be analo- 
gous to the result when the discreteness of phase space 
imposed by the Heisenberg uncertainty principle in quan- 
tum mechanics. We assume the density of states to be 
only a function of the geometrical coordination number, 
as the volume function from Eq. ([30]) reduces the degrees 
of freedom for jammed granular matter to z. The density 
of states is thereby conjectured to take the form Q: 



5(^) = (h.r 



-2d 



(33) 



where is the analogue of the Planck constant (see [8| 
for more details). 

The most populated state is the highest volume at 
z = 4 while the least populated state is the ground state 
at z = 6. This assumed form of Eq. ([33|) is an ap- 
proximation, and will be further addressed below. Since 
the term \/{hz)'^'^ is a constant, it will not influence the 
average of the observables in the partition function, al- 
though it changes the value of the entropy by a constant 
independent of (p. 

Conditions of isostaticity provide the lower bounds of 
the geometrical coordination number as Z < z, while 
considering disordered states and hard sphere conditions 
imposes ^ < 6 < 2d. which induce bounds upon the 
limits of integration in the partition function, and ac- 
count for the jamming restriction, 0jam- Here, Z is the 
mechanical coordination number related to the force bal- 
ance (that is, counting the contacts with non-zero forces) 
while z is the geometrical coordination number related to 
the geometry of the contact network. Z ranges between 
Z = d + 1, for infinite friction grains, and Z = 2d, for 
frictionless grains as given by the isostatic conjecture dis- 
cussed by Alexander [l5| and in many subsequent papers. 
Further detail on this notion is available in [a]. 

Substituting Eq. ([30]) and Eq. ([33]) and the isostatic 
condition into Eq. (|32p. we find the isostatic partition 
function: 



zux, z) = J^Ky-^^exp 



2V3' 
zX 



dz. 



(34) 



Obtaining the phase diagram is then a matter of cal- 
culating the average volume fraction, (f>{X, Z) by solving 
the partition function for different values of X and Z . 



0(X,Z) 
1 



Z;,^{X,Z) Jz Z + 2y/i 




FIG. 14: Prediction of the mesoscopic theory for (j>{X). 



The boundaries of the phase diagram are plotted in 
Fig. [2] along with the packings obtained via simulation. 
We note that Zrcp from simulations appears slightly 
higher than Z = 6, falling closer to Zrcp = 6.2 as shown 
in Fig. [21 We will use this value as Z^ax = 6.2 when con- 
sidering it in the partition function as an upper bound 
for the jamming condition in an effort to accurately ana- 
lyze the simulation results and characterize the entropy. 
Further, since the variation in Z is small as we examine 
lines of constant /z, the average of Z{p) can be used as 
the lower bound of the limit in the partition function. 



D. Compactivity 

The results presented in Fig. [71 are compared to the 
theoretical model presented above. The theoretical cal- 
culation for (t){X) is achieved exactly as described for 
(/)(Z) in Eq. ([35]) . and is presented in Fig. [TH 

In the limit of vanishing compactivity {X — > 0) for the 
theoretical model, only the minimum volume or ground 
state at z = 6 contributes to the partition function. 
Then we obtain the RCP state, (/ircp — <t>{X = Q^Z) = 
— ^—r= ~ 0.634, for all values of Z. In the limit of 

6+2\/3 ' 

infinite compactivity [X — )■ oo), the Boltzmann factor 

e~s^ 1, and the average in ([^S]) is taken over all 
the states with equal probability. Assuming I, 
the leading contribution to the average value is from the 
highest volume a,t z — Z and therefore we obtain 



4>KLv{Z) 



Z + 2V3' 



(36) 



The dotted line in Fig. [2l is a plot of the equation of 
state presented in Eq. p6p . It is an important result that 
the values of (/)rlp(-^) well fit = Z+2V3 ' other 
hand, there are states to the right of the RCP line in Fig. 
[H These states are a manifestation of the microscopic 
fluctuations, and not taken into account by the present 
mesoscopic theory. 
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The general shape of 4>{X), achieved via simulation 
well matches that as predicted by the mesoscopic the- 
ory, if not in magnitude, when using an exponential form 
for g{z). The compactivity achieves maximal value at 
the minimal available volume fraction, as observed both 
when examining the fluctuations in Voronoi volume and 
within the confines of the mesoscopic theory. Packings 
near RLP, being the least dense, have the greatest room 
to be further compacted, or increase density, and there- 
fore approach infinite compactivity. Packings near RCP, 
being the most dense, have either minimal or zero room 
for to be further compacted, and therefore cannot in- 
crease density and tend towards zero compactivity. This 
helps to establish the concept of compactivity as a static 
"effective temperature", acting as a state variable that 
may link the results of packings preparation and specific 
packing protocols to the statistical mechanics formula- 
tion. Indeed, these results for (j>{X) qualitatively resem- 
ble the compaction curves obtained in the experiments 
of ^ . We also note that the results of Fig. [T3] are qual- 
itatively similar to a simple mean-field two state model 
predicted by Edwards, where RCP and RLP are obtained 
in the limits of X — > and X oo, respectively [l3 |. 



E. Entropy 

Comparison to the theoretical model proceeds by defin- 
ing the equation of state for the entropy density. The 
entropy density is obtained as: 

{X, Z) = {w) /X + \ In Zi,o {X, Z) (37) 

This equation is obtained in analogy with equilibrium 
statistical mechanics and it is analogous to the definition 
of free energy: F = E — TS where F = —kBTlnZ is the 
free energy. We replace ksT XX, E (w). There- 
fore, F ^' E ~TS OT S ^ {E ~ F)/T = E/T + InZ is 
now s{X,Z) = {w) /X + \\nZ;so{X,Z). The partition 
function is evaluated by a numerical integration of Eq. 

for a fixed Z and as a function oi X. A numerical 
interpretation of Eq. ([55)1 then provides versus X for 
a fixed Z, a result that is plotted in Fig. [TH Using these 
two results the entropy is obtained using Eq. ([37]) . Val- 
ues of the theoretical entropy density are plotted in Fig. 
1151 for several values of Z . By fixing Z , we are equiv- 
alently imposing a fixed ^ upon the system, as Z{^) is 
determined by /i exclusively within the confines of this 
model. 

We see that the theoretical entropy density captures 
the general behavior found in the simulations as shown 
in Fig. 1131 i.e., it is maximal as we approach RLP for 
Z = A and X ^ oo while approaching the minimum en- 
tropy at RCP. Furthermore, all the curves for different 
Z approach 5 ~ InX as X ^ 0, similar to a thermal 
ideal gas. At the mesoscopic level, the entropy vanishes 
at RCP. In fact it diverges to — oo when 4> (f'RCP closer 
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FIG. 15: Prediction of the mesoscopic theory for Smoso (<?!>). 

than a constant proportional to hz (once again assum- 
ing an exponential distribution for g{z)), much like the 
Planck constant imposes a finite size in the phase space 
of quantum mechanics. Thus, we assert that the value 
of (p at which s(0) = in the theoretical model provides 
a definition of RCP at the mesoscopic level. The theo- 
retical 4i{X) arising from Fig. [14] is also in qualitative 
agreement with the simulation results of Fig. [7| 

We use hz = 0.01 in Fig. [15] such that the meso- 
scopic entropy vanishes very close to the predicted value 
of (/'RCP ~ 0.634 [8:]. However, the maximum value of 
(j)^CP « 0.642 from simulation, introduces a discrepancy 
between the theoretical and simulated models. 

The value of hz is chosen to fit the mesoscopic theory 
of Fig. [15] with simulation as close as possible, where 
the only constraint imposed by theory is hz < 1. While 
the values of both 0rcp ^-nd 0rlp(^) from the theory 
are well reproduced by the simulation results, it is clear 
that the values of s are not, as evidenced by comparing 
Fig. [T5] with Fig. [T3] For example, from simulations 
we find srlp = 2.8A aX Z = 4.3 and theory predicts 
srlp — 8.8A. This is directly due to the magnitude of hz, 
and its implications towards the density of states, g{z). 
If we examine s(X — > oo, Z), we achieve the entropy as 
a function of Z along the RLP line. When X oo, the 
equation of state in Eq. (|37p is rewritten as 

Sn,cso{X ^ c^, Z) = A In / ""^ (hzYdz, (38) 

where Z^ax = 6.2. This equation is exactly solvable, 
resulting in the following formula for the mesoscopic en- 
tropy along the RLP line. 

Wso(X^oo,Z)^Alnf - M . (39) 

Adjusting the value of hz directly affects the meso- 
scopic entropy along the RLP as defined by Eq. ([55]) . A 
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similar analysis has been used to obtain the functional 
form for (p{X ^ oo, Z) along the RLP line. 



4>{X ^ oo,Z) 



Inh, 



hz 



h? 



+ 2^3 



(hzYdz. 



(40) 



Equation (|40)) is the equation of state for the RLP line, 
and is plotted in Fig. [5] in the limit ^ when it re- 
duces to Eq. ((36l) . However, for a general hz, Eq. ([40]) 
applies. While this equation is not exactly solvable, it 
is easily seen that changing will not only impact the 
magnitude of the mesoscopic entropy along the RLP line 
of Eq. ([39|) . but also impact the values of 0rlp(2') which 
define the left most boundary of the phase diagram. Fur- 
ther, the effect on one will be the inverse of the effect on 
the other. Simply stated, within the present mesoscopic 
framework, one cannot satisfy fitting both the value of 
the entropy at, for instance, RLP and the value of (/)rlp 
in the phase diagram from simulation to the mesoscopic 
theory where hz is the only adjustable parameter. 



F. Mesoscopic and Microscopic Fluctuations 

Near RCP the mesoscopic entropy vanishes. More 
specifically, it diverges to — cx) when (/) — > 0r,cp closer 
than a constant proportional to hz, providing a char- 
acterization of RCP at the mesoscopic level. This re- 
sult qualitatively resembles behavior of the complexity 
of t he j ammed state in the replica approach to jamming 
(sol Isil l . We identify this point as a mesoscopic "Kauz- 
mann point" , in analogy with the density, or tempera- 
ture, at which the configurational entropy of a colloidal, 
or molecular, glass vanishes at the ideal glass transition 
p6| . From this point the entropy increases monotonically 
with X, being maximum for the RLP limit. An impor- 
tant result is the direct implication of a larger number of 
states available to jammed systems at RLP with respect 
to any higher volume fraction, directly implying maximal 
entropy at the RLP limit. Packings with packing frac- 
tions above RCP to (/)fcc = 0.74, the optimal packing 
fraction for spheres in 3d, do not appear in our theory 
because they exhibit some degree of order, or crystalliza- 
tion. By doing so, we explicitly do not consider crystals 
or partially crystalline packings in the ensemble. This is 
a direct consequence of setting the upper limit: z < 6. 
States with (f> > 4>ncp are new microscopic states of the 
system, and their existence requires further theoretical 
investigation. 

At RCP, we find minimal fluctuation with respect to 
Voronoi volumes associated to each grain. This implies 
the surprising conclusion of a minimal number of meso- 
scopic states for frictionless systems at RCP. Consid- 
ering these minimal fluctuations to be essentially zero, 
RCP has zero entropy and no fluctuations with respect 
to a mesoscopic coarse-graining over the ensemble. This 
is the frictionless jamming transition [J, or J-point 



[liUjQ- We see that, in principle, at the mesoscopic level 
this transition point is well-defined. However a meso- 
scopic state parameterized by a given average coordina- 
tion number contains many microscopic states which are 
averaged out in the coarse-graining procedure to calcu- 
late the volume function at the quasi-particle level. For 
instance, while the isostatic condition requires Z = 2d 
contacts per grain averaged over the entire packing, it 
makes no implication towards the exact distribution of 
contacts per individual grain. This allows for the exis- 
tence of microscopic states with grains having Z < 2d 
and Z > 2d within the packing. Therefore we expect 
that these microscopic states contribute to a nonzero en- 
tropy at the J-point at RCP. Indeed the Shannon entropy 
calculation of Section |TTlD finds the entropy of RCP to 
be srcp = l.lA. 

The full entropy should consider both mesoscopic and 
microscopic contributions, such that 



(41) 



The mesoscopic contribution is obtained via the the- 
ory, while the microscopic contribution can be obtained 
herein using the Shannon entropy method. We know that 
Srcp = 1-lA from simulations and Smcso — from the- 
ory. Therefore, the total entropy at RCP is just the mi- 
croscopic entropy, srcp = SmicroiRC P) , which is equal 
to I.IA from simulations. Thus, Smicro(^C'-P) = I.IA. 
This result is understood since all the jammed states are 
degenerate around the mesoscopic ground state with the 
coordination number z = 6. As noted, these states still 
have slightly different volume fractions, which leads to 
the microscopic fiuctuations which are coarse-grained in 
the mesoscopic theory. At the present time we do not 
have a theory to explain the value of I.IA since the meso- 
scopic theory does not include microscopic states nor fluc- 
tuation in the coordination number. Next, we make the 
additional assumption that the microscopic entropy is in- 
dependent of the volume fraction and can therefore con- 
sider SmicToi'P) = Smicro(-RC'P) = I.IA for all valucs of cj) 
between RCP and RLP. The result is the total entropy 
of the packing as 



s(0) 



,(0) + l.lA 



(42) 



for any 0, where Smcso is calculated by Eq. ([57|) . Equa- 
tion is plotted in Fig. [151 Where the mesoscopic 
entropy is augmented by its fixed value at RCP, and this 
value is srcp = 1-lA. We reiterate that in Fig. [151 the 
addition of 1.1 to Smcso{4>) makes two assumptions. First, 
its assumes that Sn-icsoiRCP) = as predicted by theory. 
Second, it assumes that Smicro(0) = Smicro(^C^') = I.IA, 
as determined by the Shannon entropy method, and is 
not an explicit function of (j). 

The mesoscopic entropy is Smcso — not only at RCP 
but for 0RCP < 4> < 4>FCC- This implies an intrinsic dif- 
ference between the current theory and Edwards' statis- 
tics given by Eq. l|77)) leading to the separation of the 
entropy in terms of the different length scales as in Eq. 
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FIG. 16: Volume fraction as a function of the inverse of com- 
pactivity, for Z = 4, as calculated using theoretical methods, 
including negative compactivity values, for ranging from 
0.001 to 0.1. When X discontinuously jumps from oo to — oo, 
4>{X) exhibits continuous behavior. However, when X con- 
tinuously goes from 0~ to 0^, 4>{X) exhibits a discontinuous 
jump from ~ 0.536 to ~ 0.634. 



G. Negative Compactivity 

Analyzing the partition function from a mathematical 
approach, we are interested in the concept of negative 
compactivity and its effect on the equation of state. 

When X ^ 0+, (t>{X 0+) 0RCp. Under the as- 
sumption of a very large g{z) at z = Z, with respect to 
any higher z, X — > +oo, (f){X — > -|-oo) — > 0rlp(^) — 
-—-^—j=^ where Z is the lower limit of integration in Eq. 

([35|) . This occurs because when using the density of 
states of Eq. ([55]) . with — > 0, g{z) ^ S{z — Z). How- 
ever, as evidenced above, an exponential form for g(z) 
may not well reproduce simulation results for the en- 
tropy. Altering g{z) is shown to better reproduce the 
entropy equation of state, but shifts the predicted value 
of 0RLP higher. Examination of a negative compactivity 
within the above presented statistical mechanics frame- 
work should allow us to achieve the minimum value of 
RLP in the limit X — > 0^. For a background on nega- 
tive temperature states in equilibrium statistical mechan- 
ics see Landau and Lifshitz book on statistical physics. 
Note that a negative temperature state is "hotter" than a 
state at absolute zero temperature and any positive tem- 
perature state. The state with T +oo is physically 
identical to the state with T —^ —oo. As in magnetic sys- 
tems where negative temperature states can be observed, 
granular matter is characterized by a bounded volume 
fraction and " thermalization" of volume at a negative 
temperature is possible, in principle. Next, we analyze 
those states with a negative compactivity. 

When AT ^ 0+, the Boltzmann factor in the parti- 
tion function of Eq. ([M)) tends towards zero. As such, 
the largest value of z, or the smallest value of 1/z will 
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FIG. 17: Entropy plotted for Z = 4, as calculated using theo- 
retical methods, including negative compactivity values. S{<l)) 
tends towards — oo at <^ = 0.536 {X 0~) and (j) = 0.634 
{X ^0^), the minimal value of RLP and RCP, respectively, 
for any value of h^. Dashed vertical lines show 4i{X) for 
X ±cx), acknowledging the lowest physically achievable 
volume fraction ai Z = 4 for a particular value of hz. Larger 
values of hz result in larger values of ±co). Following 

[20| we call (f){X 0~) — » 0vlrp (very loose random pack- 
ing), (f){X ±oo) 0RLP and (j){X O"*") ?!>rcp. In the 
limit of /iz — > the difference between VLRP and RLP van- 
ishes =J> '^'^'^^^o ~^ fl^RLP. The maximum entropy is always at 
RLP. Negative compactivity states are very difficult to obtain 
with current protocols. 

give the largest value of the Boltzmann factor when us- 
ing the partition function to calculate observable aver- 
ages. This results in the calculation of the RCP state. 
However, when X —^ 0~ , the Boltzmann factor in the 
partition function tends towards infinity, not zero, and 
the largest value of Z, or the smallest value of 1/z will 
give the largest value of the Boltzmann factor. When cal- 
culating the average volume function in either case, the 
density of states will not greatly impact the results with 
respect to the contribution from the Boltzmann factor. 
Therefore, when X —^ 0~ , the average volume fraction 
will reduce to the predicted value of the RLP line where 
^ 0, </)(A ^ 0-,Z) = 0RLp(^) = Figure [H 

exemplifies this phenomenon, plotting (j){X) aX Z ~ A for 
different values of hz- 

The entropy equation of state should achieve the same 
values as A — > +00 and A — > — 00, since the Boltzmann 
factor approaches unity in either case. Although not as 
obvious, the same can be said for A ^ 0^ and A — > 0^. 
The equation of state (|37p can be rewritten as: 

Wso(A, Z)/\ = In / V,)^-2'*e- At(-(-)-<-(-)»dz. 
J z 

(43) 

In the case of A ^ 0+, A is positive, and w{z) > 
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{w{z)), as RCP represents the lowest attainable value 
of {w{z)). The exponentiated term in the partition func- 
tion is always negative, such that the entropy approaches 
—oo. Conversely, when X —f 0~ , X is negative, and 
w{z) < (wlz)), and we obtain the highest attainable 
value of Again, the exponentiated term in the 

partition function is always negative, such that the en- 
tropy approaches — oo. 

Figure [T7| displays the equation of state for Z = 4 as 
calculated by the theoretical model of Eq. (j43|) . We 
use Z = A and Zmax = 6 as the limits of integration 
for the partition function of Eq. . The value of hz is 
increased to show that larger values of 5(2) at RLP result 
in a more precise agreement between (j){X ±00) and 
(f){X 0^ ) . The dashed vertical lines of Fig. [T7| show 
for X ±00, acknowledging the lowest physically 
achievable volume fraction at Z = A for a particular value 
of hz- Larger values of hz result in larger values of (j){X 
±00), allowing for the existence of a greater range of 
packings with a negative compactivity, having (j){X — s- 
0") < < 0(X ^ -00). 

Some models exist where the concept of a negative tem- 
perature finds physical meaning, including nuclear spins 
and semiconducting lasers. In Ref. [20j an attempt is 
made to include the concept of a negative temperature 
within a statistical mechanics framework. By means of a 
lattice model in 2d, utilizing a discrete phase space, me- 
chanically stable packings, or microstates, are shown to 
exist beneath the volume fraction with the largest num- 
ber of microstates, at a particular fi. The highest entropy 
occurs when the largest number of microstates are avail- 
able, and is the equivalent of (/jrlp in the above presented 
mesoscopic theory. Under the assumption of an ergodic 
exploration of the volume fractions available to the lattice 
model, this implies that packings with (j) < 0rlp exist, 
with entropy below the maximal value, and that these 
packings can be explained via the concept of a negative 
temperature. Reference f20| thereby considers (?!)rlp to 
be the " loosest possible random packing that is mechan- 
ically stable that one can achieve by pouring grains". 
Below this limit there exists RVLP "random very loose 
packings" with negative temperature. 

Although the present work analyzes packings in 3d, 
the results of mesoscopic theory find agreement with the 
simulations of [l^j. The main results of Q, Fig. [2] and for 
instance Eq. (|36p have been obtained in the limit of hz — > 
0. If we relax this constraint then negative compactiv- 
ity states are possible. We assert that these states exist 
because there is an upper bound in the volume function 
at z = Z. Systems with unbounded volume functions 
(Hamiltonians) do not allow for negative compactivity 
(temperature) states. In this context we have the def- 
inition of the following limits according to the entropy 
and compactivity: The RCP limit is (j){X ^ 0+) = 0rcp 
and minimum entropy: Neglecting the negative diverg- 
ing of the entropy we have S'rcp — S{X 0+) 0. 
The RLP limit is defined as the maximum entropy in 
the limit (j){X — > ±00) = <^r,lp- If hz is finite then 



VLRP appears as (f>{X ^ ) = (^vlrp and minimum 
entropy. Again, neglecting the divergency, we obtain 
S-VLRP = S{X ^ 0-) ^ 0. In the limit of hz ^ 
the difference between RLP and VLRP vanishes and we 
have only one weU defined RLP as: ^ 0vlrp- 

Fig. 1171 shows a maximal entropy at 0rlp, indicat- 
ing the largest number of available microstates to the 
system. The introduction of a negative compactivity, as 
described above, allows the theory to probe states such 
that (f) < (t>B.LP- It becomes apparent that the range of 
(j) in which these states may exist is directly related to 
the magnitude of hz, decreasing as the discretization of 
phase space within the confines of the mesoscopic theory 
such that in the limit of a continuous phase space none of 
these packings are mechanically stable. Our theoretical 
model includes the concept of a discrete phase space for 
jammed grains, and our simulations show that hz < 1, 
but not necessarily /i^ <C 1. Further, the "split" algo- 
rithm utilized simulates a pouring of grains with respect 
to the method of packing creation. It is possible that the 
concept of a negative compactivity could help to explain 
areas in the phase diagram of Fig. [5] unavailable within 
the scope of the present study. 



IV. OUTLOOK 

Although extensive detail is presented in this study re- 
garding the various steps necessary to analyze the equa- 
tions of state, several questions remain unclear. 

(a) The derivation of the entropy for jammed granular 
matter is explicitly calculated via fluctuations in Voronoi 
cell volume for each packing presented herein. As we are 
studying random packings, crystal states are not achieved 
with anything other than measure zero probability. Thus, 
the highest available volume fraction for any given pack- 
ing is RCP, 4> ~ 0.64, not FCC, (j> = 0.74, as corroborated 
by simulation results. This result is the apparent limit 
of the preparation protocols used herein. It can be said 
with certainty that a FCC packing has zero fluctuation 
with respect to their constituent Voronoi cells. At RCP 
simulations reveal a non-zero fluctuation, as discussed 
above. This approach, however, does not account for 
the possibility that for packings above RCP, but below 
FCC, may achieve a continuous, monotonic, degree of 
crystallization. Such a condition would permit a contin- 
uous decrease in fluctuation to exactly zero, from RCP 
to FCC, which can be studied within the scope of a more 
complete theory that includes random and partially crys- 
tallized packings. 

Below we elaborate on a possible scenario to rationalize 
the transition from disorder at RCP to order at FCC. 
If the microscopic fluctuations are not subtracted from 
RCP then the compactivity curves presented in Fig. [7| 
no longer reach a plateau when approaching RCP, but 
reach a finite, non-zero, value. This opens the possibility 
that a true thermodynamic phase transition may occur at 
RCP between a disordered phase and an ordered phase. 
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It remains possible that a phase transition occurs at RCP, 
and packings of higher volume fraction need not preserve 
the properties of a fully random system. It remains an 
open topic how one would define compactivity between 
RCP and FCC, as compactivity as been herein attributed 
to packing protocols resulting in random packings. 

The mesoscopic theory of 8] utilized herein considers 
a vanishing entropy at, or near, RCP. Taking into consid- 
eration the FCC state, a more complete theory including 
packings between RCP and FCC should be character- 
ized, such that the entropy of jammed matter approaches 
zero when approaching FCC, not RCP. Furthermore, the 
result that X ^ at RCP is merely an artifact of the 
mesoscopic theory that neglects microscopic fluctuations. 
A full theory would obtain X ^ at FCC. Figure [H] dis- 
plays a possible interpretation for an extension of the en- 
tropic equation of state. The entropy attains a maximal 
value at 0rlp, as predicted by the existing mesoscopic 
theory. When 0rlp < <t> < 0rcp i the packing consists of 
purely random states, and the entropy decreases as we 
approach RCP from lower volume fractions, also as pre- 
dicted. At some point close to <?!>rcp, the entropy deviates 
from its predicted decrease to zero at (/)rlp, and follows 
a different branch. When (/)rcp < (j) < 4>fcc a coexis- 
tence between random and crystalized microstates may 
exist, ultimately leading to a purely crystalized packing 
at FCC. The exact nature of the transition, continuous 
or discontinuous, from purely random states to a coexis- 
tence of states remains an open topic. The incorporation 
of microscopic fluctuations and microscopic crystalized 
states into the existing mesoscopic theory may result in a 
more complete characterization of the entropy of jammed 
granular matter, a work currently in progress. 

(b) An additional assumption made in calculating the 
compactivity of jammed granular matter is that those 
packings along the RLP line have X = oo. This as- 
sumption directly allows for the (/>(X) equation of state 
without any constants due to integration techniques. It 
appears to be reasonable to presume that X is large along 
the simulated RLP line, with respect to X of all other 
packings used herein, and that the constant term -^rrr — v 
is very small. 

(c) As mentioned above, there exists a small fluctua- 
tion density in Voronoi volume at RCP. The fluctuation 
density increases as we increase n, and comes to a max- 
imal value at a particular range of n, where n is larger 
as we approach lower values of (jj as shown in Fig. [5l 
This increase in value is due to selecting of clusters of 
n Voronoi volumes, as opposed to n randomly chosen 
grains. In the case of randomly chosen grains, the fluctu- 
ation density remains constant, equal to the fluctuation 
at n = 1. Such a change implies a correlation between 
the Voronoi volumes in clusters of grains, exemplified by 
Eq. ([7]). This correlation may create a scale separation, 
such that microscopic correlations dominate the fluctua- 
tions at lower n and mesoscopic correlations dominate at 
higher n. This scale separation would result in a differ- 
ence between local and global compactivities, suggested 




FIG. 18: A possible extension of the entropy equation of state 
to include packings between RCP and FCC. For (^rlp < 4> < 
4'B.CP, the packing consists of only random states, and the 
entropy decreases when approaching RCP from lower volume 
fractions, as predicted by the mesoscopic curve. At <jf>RCP, 
the entropy does not decrease to zero, as predicted by the 
mesoscopic curve, and follows a different branch, achieving 
exactly zero entropy at </>fcc- When <^r,cp < (f> < 4>fcc a 
coexistence between random and crystalized microstates may 
exist. The possible existence of a transition that may define 
RCP remains an open question. 



in [47[ , and is a topic of continuing study. It also remains 
possible that fluctuation densities calculated using larger 
n sized clusters overuse the data set by largely repeating 
constituent Voronoi cells in cluster volumes, thereby ren- 
dering the fluctuation densities questionable for packings 
approaching RLP. Again, it should be noted that simply 
using the microscopic fluctuations does not greatly effect 
the entropy calculations. 

(d) The mesoscopic theory predicts S — > — oo as we 
approach RCP for any /i. This concept does not adhere 
to physical measurements, where the minimal entropy 
should be zero. Taking this into consideration, another 
method must be used to calculate the entropy of packings 
at RCP, independent of the distribution of Voronoi vol- 
umes used to facilitate the mesoscopic theory. The Shan- 
non entropy calculation applies graph theoretical meth- 
ods, resulting in a non-zero value for the entropy of RCP, 
well suited for a more complete equation of state. While 
it can be said with some degree of certainty that the 
Shannon entropy calculation contributes to the entropy 
at RCP, it is unclear whether it does so completely. Other 
methods may be available that provide all of the entropy 
at RCP, or give additional terms to the Shannon entropy, 
another topic of continuing study. Further, as discussed 
above with respect to fluctuation density, packings be- 
tween RCP and FCC have partial degrees of crystalliza- 
tion, such that the entropy of an FCC is ex actly zero, or 
— oo as described in the theory. Recent work [23| suggests 
that the entropy experiences an increase immediately fol- 
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lowing RCP, due to mixing between random and crystal 
states. Incorporation of these ideas into the present work 
remains a topic of ongoing study. 

(e) While the concept of negative compactivity works 
well as a mathematical tool for achieved the VLRP state 
of jammed granular matter, it remains difficult to at- 
tribute a physical meaning to such a condition. Recent 
studies in 20] suggest that negative compactivity probes 
mechanically stable states that exist beneath RLP, but 
are not accessible by means of grain pouring. It should 
be noted that X < states may not be thermalized with 
X > 0. Such ideas will remain the topic of future studies. 

(f) An exponential density of states with a very small 
value of hz provides a very accurate phase diagram, in 
comparison with simulations, but values for the entropy 
that are different by a factor of 4 times larger than those 
achieved via simulation. It remains possible that the 
methods presented herein simply do not capture all of the 
entropy for jammed matter, and larger values from sim- 
ulation are possible, such that corrections to the meso- 
scopic theory could be minimized. Further, a new the- 
ory that carefully analyzes all packings up from 0r,lp 
through 0FCC =0.74 may result in the ability to capture 
all of the entropy for a given packing. Current work is 
approaching this idea from a microscopic level. 

(g) The differences displayed between theoretical and 
simulated entropy can to some extent be considered 
within the scope of disagreement between classical and 
quantum entropy. Classical entropy measurements are 
primarily interested in AS", having a minimum value 
of when approaching the ground state, though S 

— oo. Quantum entropy measurements assume a mini- 
mum value of phase space over which one can integrate 
degrees of freedom for a give system. This results in a 
S = exactly at the ground state. Determination of 
a well defined minimum phase space volume for jammed 
matter would adjust the number of available microstates, 
f2, within the micro-canonical ensemble as presented by 
Edwards. Such a model may be available through a care- 
ful analysis of the microscopic entropy, including par- 
tially crystalized states above c/jkcp through the crystal 
state of ^Fcc , the highest achievable packing fraction for 
identical spheres in 3d. 

V. CONCLUSIONS 

Simulation results as derived from mesoscopic fluctu- 
ations are presented in Figs. [7] and [51 When comparing 
all the packings with different Z{^) and (f>, the maximum 
entropy is at the minimum volume fraction of RLP (f>m^p 
when X ^ oo but only infinite friction. The minimum 
entropy is found for the RCP state at (pucp for X ^ 0, 
now for all the values of friction, indicating the degener- 



acy of the RCP state. It is commonly believed that the 
RCP limit corresponds to a state with the highest number 
of configurations and therefore the highest entropy. This 
belief is expressed for instance in the definition of RCP 
as the maximally random jammed state (l^ . However, 
here we show that the states with a higher compactivity 
have a higher entropy, corresponding to the looser RLP 
packings. Within a statistical mechanics framework of 
jammed matter, this result is a natural consequence, and 
gives support to such an underlying statistical picture. 

Each curve in the Figs. [T3] and [15] correspond to a 
system with a different Z^jj), as calculated using a meso- 
scopic ensemble. When comparing all the packings, the 
maximum entropy is at 0rlp and X ^ oo while the en- 
tropy is minimum for 0rcp at X — > 0''". Following the 
Z = 4 line in the phase diagram we obtain the entropy for 
infinitely rough spheres showing a larger entropy for the 
RLP than the RCP. The same conclusion is obtained for 
the other packings at finite friction (4 < Z{y) < 6). We 
conclude that the RLP states are more disordered than 
the RCP states. Approaching the frictionless J-point, 
/i — > (Z = 6) the entropy vanishes. More precisely, 
it vanishes for a slightly smaller (f> than (/ircp of the or- 
der hz- Strictly speaking it diverges to — oo at 0rcp as 
S InX for any value of Z, in analogy with the classical 
equation of state. However, this is an unphysical limit, as 
it would be considering distances in phase space smaller 
than the minimal distance in the jamming phase space. 
Thus we consider only packings with an entropy density 
greater than or equal to as "physical" packings. We 
note that the compactivity curves from the theoretical 
model match simulation with accuracy. The theoretical 
entropy fails to agree with the entropy from simulation 
in magnitude, but reproduces the overall shape. While 
increasing hz such that the magnitude of the entropy 
from the mesoscopic theory would decrease, the RLP line 
would no longer be well reproduced. As simulations and 
theory are in strong agreement with respect to the RLP 
line, increasing hz does not appear to be a reasonable 
amendment to the mesoscopic partition function. 

In summary, a notion of disorder is presented that ap- 
plies to frictional hard spheres, as well as frictionless ones. 
The entropy reveals interesting features of the RCP and 
RLP states such as the fact that RLP is maximally ran- 
dom with respect to RCP and that both limits can be 
defined in terms of the entropy and equation of state. 
Overall, the agreement between theory and simulation is 
sufficient to indicate that the methods presented herein 
are appropriate for evaluating the entropy of jammed 
matter. 
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